/*
 * AnalyzeCharges.h
 *
 *  Created on: May 26, 2011
 *      Author: marchi
 */

#ifndef ANALYZECHARGES_H_
#define ANALYZECHARGES_H_

#include "Metric.h"
#include "RhoAverage.h"
#include "RhoDipAverage.h"
#include <cstdlib>
#include <cstdio>
#include "time.h"
#include "pdbio.h"
#include "confio.h"
#include "CellSymmetry.h"
#include "EpsofK.h"
#include "Grid_SpecTemplates.hpp"

//enum ElecType {Slt, Sol, Tot};

static Rho * ro=NULL;
static RhoDip * rodp0=NULL;
static int nframes=0;
static int mycounter=0;
static CellSymmetry * symmetry=NULL;
static FILE * fg3, * fg4;

static int Analyze(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
    /* Here, you can do whatever analysis your program requires for a frame. */
	Averages * Avgs=static_cast<Averages *> (data);

	RhoAverage ** roavg0=Avgs->Ro;
    RhoAverage * & roavg=*roavg0;
    RhoDipAverage * rodipavg0=Avgs->RoDip;
    RhoDipAverage & rodipavg=*rodipavg0;
    if(mycounter == 0) {
    	fg3=fopen("MyEpsilon_in.xvg","w+");
    	roavg[Sol][0]=0.0;
    }
    mycounter++;
    srand (time(NULL) );

    Grid<1>::set(fr->box);
	Grid<DIM>::set(fr->box);


    if(!symmetry) symmetry=new CellSymmetry(top,fr,pbc);
    if(!ro) {
    	ro=new Rho [ngrps];
    }

    Metric Met;
    Met(fr->box);
    symmetry->Center();

    for(int i=0;i<ngrps;i++)
    	chg[i].setCoord(Met,fr->x,cidx[i]);

    MyDip.setCoord(Met,fr->x,cidx[Sol]);


    rodipavg.Density(MyDip);

    roavg[Slt].Density(chg[Slt]);
    roavg[Sol].Density(chg[Sol]);


    nframes++;
    return 0;
}
static int AnalyzeFluctuations(t_topology *top, t_trxframe *fr, t_pbc *pbc,
              int nr, gmx_ana_selection_t *sel[], void *data)
{
    /* Here, you can do whatever analysis your program requires for a frame. */
	Averages * Avgs=static_cast<Averages *> (data);

    RhoAverage * sofk0=Avgs->RhoSqrK;
    RhoAverage & sofk=*sofk0;
    RhoAverage * sofkdip0=Avgs->RhoSqrK_Dip;
    RhoAverage & sofkdip=*sofkdip0;

    if(mycounter == 0) {
    	fg3=fopen("MyEpsilon_in.xvg","w+");
    	fg4=fopen("MyLinear_in.xvg","w+");
    }
    mycounter++;
    srand (time(NULL) );

    Grid<1>::set(fr->box);
	Grid<DIM>::set(fr->box);


    if(!symmetry) symmetry=new CellSymmetry(top,fr,pbc);
    if(!ro) {
    	ro=new Rho [ngrps];
    }
    if(!rodp0) {
    	rodp0=new RhoDip;
    }
    Metric Met;
    Met(fr->box);
    symmetry->Center();
    for(int i=0;i<ngrps;i++)
    	chg[i].setCoord(Met,fr->x,cidx[i]);
    MyDip.setCoord(Met,fr->x,cidx[Sol]);

    RhoDip & rodp=*rodp0;

    rodp.Density(MyDip);
    ro[Sol].Density(chg[Sol]);
    ro[Slt].Density(chg[Slt]);
    ro[Tot]=ro[Sol]+ro[Slt];

    EpsofK MyEps_rho,MyEps_dip;

    MyEps_rho(ro[Sol]);
    MyEps_dip(rodp);
    sofk+=MyEps_rho;
    sofkdip+=MyEps_dip;
    sofk.Accumulate();
    sofkdip.Accumulate();
    if(nframes && !(nframes%100)) {
    	Rho MyEpsRho;
    	Rho MyEpsDip;
    	Rho rosol;
    	MyEpsRho=sofk.getavg();
    	MyEpsRho.RdfK(fg3,"S(k) by Rho(k)");

    	rosol.Linear(MyEpsRho,ro[Slt]);
    	double Myx[DIM]={Grid<1>::getCO()[0][0]*0.5,Grid<1>::getCO()[1][1]*0.5,Grid<1>::getCO()[2][2]*0.5 };

    	rosol.Rdf(fg4,Myx,"rho_sol (r) Linear");

//    	MyEpsDip=sofkdip.getavg();
//    	MyEpsDip.RdfK(fg3,"S(k) by M(k)");
    }

    nframes++;
    return 0;
}

#endif /* ANALYZECHARGES_H_ */
